Statistical-physics-based reconstruction in compressed sensing 



F. Krzakala \ M. Mezard ^ F. Sausset 2, Y. F. Sun^-^ and L. Zdeborova ^ 
^ CNRS and ESP CI PansTech, 10 rue Vauquehn, UMR 7083 Gulliver, Pans 75005, France. 
2 Umv. Paris-Sud & CNRS, LPTMS, UMR8626, Bat. 100, 91405 Orsay, France. 
^ LMIB and School of Mathematics and Systems Science, Beihang University, 100191 Beijing, China. 
* Institut de Physique Theorique, IPhT, CEA Saclay, 
and URA 2306, CNRS, 91191 Gif-sur-Yvette, France^ 
(Dated: June 7, 2012) 

Compressed sensing is triggering a major evolution in signal acquisition. It consists in sampling 
a sparse signal at low rate and later using computational power for its exact reconstruction, so that 
only the necessary information is measured. Currently used reconstruction techniques are, however, 
limited to acquisition rates larger than the true density of the signal. We design a new procedure 
which is able to reconstruct exactly the signal with a number of measurements that approaches 
the theoretical limit in the limit of large systems. It is based on the joint use of three essential 
ingredients: a probabilistic approach to signal reconstruction, a message-passing algorithm adapted 
from belief propagation, and a careful design of the measurement matrix inspired from the theory 
of crystal nucleation. The performance of this new algorithm is analyzed by statistical physics 
methods. The obtained improvement is confirmed by numerical studies of several cases. 

The ability to recover high dimensional signals using only a limited number of measurements is crucial in many 
fields, ranging from image processing to astronomy or systems biology. Example of direct applications include speeding 
up magnetic resonance imaging without the loss of resolution, sensing and compressing data simultaneously [T and 
the single-pixel camera [2] ■ Compressed sensing is designed to directly acquire only the necessary information about 
the signal. This is possible when the signal is sparse in some known basis. In a second step one uses computational 
power to reconstruct the signal exactly [11 [31 S] . Currently, the best known generic method for exact reconstruction is 
based on converting the reconstruction problem into a convex optimization one, which can be solved efficiently using 
linear programming techniques |3l [4j . The £1 reconstruction is able to reconstruct accurately, provided the system 
size is large and the ratio of the number of measurements M to the number of non-zeros K exceeds a specie limit 
which can be proven by careful analysis [HIS]. However, the limiting ratio is signicantly larger than 1. In this paper 
we improve on the performance of £1 minimization, and in the best possible way: we introduce a new procedure that 
is able to reach the optimal limit M/K — >■ 1. This procedure, which we call seeded Belief Propagation (s-BP) is based 
on a new, carefully designed, measurement matrix. It is very powerful, as illustrated in Fig. [T] Its performance will 
be studied here with a joint use of numerical and analytic studies using methods from statistical physics [H]. 

Reconstruction in Compressed Sensing 

The mathematical problem posed in compressed-sensing reconstruction is easily stated. Given an unknown signal 
which is a iV-dimensional vector s, we make M measurements, where each measurement amounts to a projection of s 
on some known vector. The measurements are grouped into a A/-component vector y, which is obtained from s by a 
linear transformation y = Fs. Depending on the application, this linear transformation can be for instance associated 
with measurements of Fourier modes or wavelet coefficients. The observer knows the M x N matrix F and the M 
measurements y, with M < N. His aim is to reconstruct s. This is impossible in general, but compressed sensing 
deals with the case where the signal s is sparse, in the sense that only K < N oi its components are non-zero. We 
shall study the case where the non-zero components are real numbers and the measurements are linearly independent. 
In this case, exact signal reconstruction is possible in principle whenever M > K + 1, using an exhaustive enumeration 
method which tries to solve y = Fx for all (^) possible choices of locations of non-zero components of x: only one 
such choice gives a consistent linear system, which can then be inverted. However, one is typically interested in large 
instances where TV ^ 1, with M = aN and K = poN. The enumeration method solves the compressed sensing 
problem in the regime where measurement rates are at least as large as the signal density, a > po, but in a time which 
grows exponentially with N, making it totally impractical. Therefore a — pq is the fundamental reconstruction limit 
for perfect reconstruction in the noiseless case, when the non-zero components of the signal are real numbers drawn 
from a continuous distribution. A general and detailed discussion of information-theoretically optimal reconstruction 
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FIG. 1: Two illustrative examples of compressed sensing in image processing. Top: the original image, the Shepp-Logan 
phantom, of size A*' — 128^, is transformed via one step of Haar wavelets into a signal of density po ~ 0.15. With compressed 
sensing one is thus in principle able to reconstruct exactly the image with M > poN measurements, but practical reconstruction 
algorithms generally need M larger than poN. The five columns show the reconstructed figure obtained from M = aN 
measurements, with decreasing acquisition rate a. The first row is obtained with the ^i-reconstruction algorithm [3] [4 for a 
measurement matrix with iid elements of zero mean and variance The second row is obtained with belief propagation, using 

exactly the same measurements as used in the £i reconstruction. The third row is the result of the seeded belief propagation, 
introduced in this work, which uses a measurement matrix based on a chain of coupled blocks, see Fig.|4] The running time of 
all the three algorithms is comparable (asymptotically they are all quadratic in the size of the signal). In the second part of 
the figure we took as the sparse signal the relevant coefficients after two-step Haar transform of the picture of Lena. Again for 
this signal, with density po = 0.24, the s-BP procedure reconstructs exactly down to very low measurement rates (details are 
in Appendix G, the data is available online [7]). 



has been developed recently in [SUIO]. 

In order to design practical 'low-complexity' reconstruction algorithms, it has been proposed [31 13] to find a vector 
X which has the smallest ii norm, X^iLi within the subspace of vectors which satisfy the constraints y = Fx, 
using efficient linear programming techniques. In order to measure the performance of this strategy one can focus 
on the measurement matrix F generated randomly, with independent Gaussian-distributed matrix elements of mean 
zero and variance and the signal vector s having density < po < 1- The analytic study of the £i reconstruction 
in the thermodynamic limit N ^ oo can be done using either geometric or probabilistic methods O lllj , or with the 
replica method [T2HI1]- It shows the existence of a sharp phase transition at a value a£^{po). When a is larger than 
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this value, the £i reconstruction gives the exact result x = s with probability going to one in the large N limit, when 
a < (po) the probability that it gives the exact result goes to zero. As shown in Fig. [2j (po) > po and therefore 
the £i reconstruction is suboptimal: it requires more measurements than would be absolutely necessary, in the sense 
that, if one were willing to do brute-force combinatorial optimization, no more than pqN measurements are necessary. 

We introduce a new measurement and reconstruction approach, s-BP, that allows to reconstruct the signal by a 
practical method, which needs only w poN measurements. We shall now discuss its three ingredients: 1) a probabilistic 
approach to signal reconstruction, 2) a message-passing algorithm adapted from belief propagation [15', which is a 
procedure known to be efficient in various hard computational problems [161 117j , and 3) an innovative design of the 
measurement matrix inspired from the theory of crystal nucleation in statistical physics and from recent developments 
in coding theory [T8U21j . Some previous works on compressed sensing have used these ingredients separately. In 
particular, adaptations of belief propagation have been developed for the compressed sensing reconstruction, both in 
the context of £i reconstruction [TI] [551 US] , and in a probabilistic approach [21] . The idea of seeding matrices in 
compressed sensing was introduced in |25j . It is however only the combined use of these three ingredients that allows 
us to reach the a — po limit. 




FIG. 2: Phase diagrams for compressed sensing reconstruction for two different signal distributions. On the left-hand side 
the poN non-zero components of the signal are independent Gaussian random variables with zero mean and unit variance. On 
the right-hand side they are independent ±1 variables. The measurement rate is a = M/N. On both sides we show, from top 
to bottom: (a) The phase transition ae^ for £i reconstruction [5l 1111 [12] (which does not depend on the signal distribution) . 

(b) The phase transition q:em-bp for EM-BP reconstruction based for both sides on the probabilistic model with Gaussian (j). 

(c) The data points which are numerical reconstruction thresholds obtained with the s-BP procedure with L — 20. The point 
gives the value of a where exact reconstruction was obtained in 50% of the tested samples, the top of the error bar corresponds 
to a success rate of 90%, the bottom of the bar to a success of 10%. The shrinking of the error bar with increasing A'^ gives 
numerical support to the existence of the phase transition that we have studied analytically. These empirical reconstruction 
thresholds of s-BP are quite close to the a = po optimal line, and get closer to it when increasing A'^. The parameters used in 
these numerical experiments are detailed in Appendix E. (d) The line a = po that is the theoretical reconstruction limit for 
signals with continuous <^o- An alternative presentation of the same data using the convention of Donoho and Tanner _5_ is 
shown in Appendix F. 



A probabilistic approach 

For the purpose of our analysis, we consider the case where the signal s has independent identically distributed (iid) 
components: Pq{s) = ni=i[(l ~ Po)S{si) + po4>o{si)], with < po < 1. In the large-iV limit the number of non-zero 
components is poN. Our approach handles general distributions 0o(si)- 

Instead of using a minimization procedure, we shall adopt a probabilistic approach. We introduce a probability 
measure P(x) over vectors x e which is the restriction of the Gauss-Bernoulli measure P(x) = Y[iLi{i^~ P)^i^i)~^ 
p(f>{xi)] to the subspace |y — Fx| — 26J. In this paper, we use a distribution (pix) which is a Gaussian with mean x 
and variance ct^, but other choices for (l){x) are possible. It is crucial to note that we do not require a priori knowledge 
of the statistical properties of the signal: we use a value of p not necessarily equal to pQ, and the </> that we use is not 
necessarily equal to (po- The important point is to use p < 1 (which reflects the fact that one searches a sparse signal). 
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Assuming that F is a random matrix, either where all the elements are drawn as independent Gaussian random 
variables with zero mean and the same variance, or of the carefully-designed type of 'seeding matrices' described 
below, we demonstrate in Appendix A that, for any po-dense original signal s, and any a > po the probability P{s) 
of the original signal goes to one when N —> oo. This result holds independently of the distribution 0o of the original 
signal, which does not need to be known. In practice, we see that s also dominates the measure when N is not very 
large. In principle, sampling configurations x proportionally to the restricted Gauss-Bernoulli measure P(x) thus gives 
asymptotically the exact reconstruction in the whole region a > pq. This idea stands at the roots of our approach, 
and is at the origin of the connection with statistical physics (where one samples with the Boltzmann measure) . 
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FIG. 3: When sampling from the probability f'(x), in the limit of large A'^, the probability that the reconstructed signal x is at 
a squared distance D = ~ Si)^ /N from the original signal s is written as ce^*'^' where c is a constant and 9{D) is the 

free entropy. Left: $(-D) for a Gauss-Bernoulli signal with po = p — 0.4 in the case of an unstructured measurement matrix F 
with independent random Gaussian-distributed elements. The evolution of the EM-BP algorithm is basically a steepest ascent 
in $(-D) starting from large values of D. It goes to the correct maximum at D = for large value of a but is blocked in 
the local maximum that appears for a < QfEM-Bp(po ~ 0.4) « 0.59. For a < po, the maximum is not at D = and exact 
inference is impossible. The seeding matrix F, leading to the s-BP algorithm, succeeds in eliminating this local maximum. 
Right: Convergence time of the EM-BP and s-BP algorithms obtained through the replica analysis for po — 0.4. The EM-BP 
convergence time diverges as a — >■ oem-bp with the standard L = 1 matrices. The s-BP strategy allows to go beyond the 
threshold: using ai — 0.7 and increasing the structure of the seeding matrix (here L — 2,5,10,20), we approach the limit 
a = Po (details of the parameters are given in Appendix E) . 



Sampling with Expectation Maximization Belief Propagation 



The exact sampling from a distribution such as P(x), eq. (A2|, is known to be computationally intractable |27j . 
However, an efficient approximate sampling can be performed using a message-passing procedure that we now describe 
[^HdHHSS]. We start from the general behef-propagation formalism [T71 [3TJ [32] : for each measurement p = I, . . . ,M 
and each signal component i = 1,...,N, one introduces a 'message' mi^^{xi) which is the probability of Xi in 
a modified measure where measurement p has been erased. In the present case, the canonical belief propagation 
equations relating these messages can be simplified [HI [22H241 [33] into a closed form that uses only the expectation 
af\^ and the variance vf2,^ of the distribution mf\^{xi) (see Appendix B). An important ingredient that we add 
to this approach is the learning of the parameters in P(x): the density p, and the mean x and variance ct^ of the 
Gaussian distribution (/'(x). These are three parameters to be learned using update equations based on the gradient 
of the so-called Bethe free entropy, in a way analogous to the expectation maximization This leads to the 

Expectation Maximization Belief Propagation (EM-BP) algorithm that we will use in the following for reconstruction 
in compressed sensing. It consists in iterating the messages and the three parameters, starting from random messages 
o-i^fj^ and v^^^, until a fixed point is obtained. Perfect reconstruction is found when the messages converge to the 
fixed point ai^u = Si and Vi^^ = 0. 
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Like the £i reconstruction, the EM-BP reconstruction also has a phase transition. Perfect reconstruction is achieved 
with probabihty one in the large- TV limit if and only if a > aEM-BP- Using the asymptotic replica analysis, as 
explained below, we have computed the line Q:em-bp(/Oo) when the elements of the M x N measurement matrix F are 
independent Gaussian random variables with zero mean and variance 1/N and the signal components are iid. The 
location of this transition line does depend on the signal distribution, see Fig. [2j contrary to the location of the £i 
phase transition. 

Notice that our analysis is fully based on the case when the probabilistic model has a Gaussian 0. Not surprisingly, 
EM-BP performs better when </> = (/)o, see the left-hand side of Fig. [2] where EM-BP provides a sizable improvement 
over ii. In contrast, the right-hand side of Fig. [2] shows an adversary case when we use a Gaussian to reconstruct 
a binary (pQ, in this case there is nearly no improvement over £i reconstruction. 

Designing seeding matrices 

In order for the EM-BP message-passing algorithm to be able to reconstruct the signal down to the theoretically 
optimal number of measurements a = po, one needs to use a special family of measurement matrices F that we 
call 'seeding matrices'. If one uses an unstructured F, for instance a matrix with independent Gaussian-distributed 
random elements, EM-BP samples correctly at large a, but at small enough a a metastable state appears in the 
measure P{^), and the EM-BP algorithm is trapped in this state, and is therefore unable to find the original signal 
(see Fig. |3|, just as a supercooled liquid gets trapped in a glassy state instead of crystallizing. It is well known in 
crystallization theory that the crucial step is to nucleate a large enough seed of crystal. This is the purpose of the 
following design of F. 

We divide the N variables into L groups of N/ L variables, and the M measurements into L groups. The number 
of measurements in the p-th group is Mp = UpN/L, so that M = X]p=i "^p] N = a N. We then choose the 

matrix elements independently, in such a way that, if i belongs to group p and /i to group q then is a random 
number chosen from the normal distribution with mean zero and variance Jq,p/N (see Fig. |4|. The matrix Jg p is 
a, L X L coupling matrix (and the standard compressed sensing matrices are obtained using L = 1 and ai ~ a). 
Using these new matrices, one can shift the BP phase transition very close to the theoretical limit. In order to get an 
efficient reconstruction with message passing, one should use a large enough ai. With a good choice of the coupling 
matrix Jp,q, the reconstruction first takes place in the first block, and propagates as a wave in the following blocks 
p — 2,3, ... , even if their measurement rate Up is small. In practice, we use a2 = ■■■ = ct^ ^ a' , so that the total 
measurement rate is a = [ai + {L— l)a']/L. The whole reconstruction process is then analogous to crystal nucleation, 
where a crystal is growing from its seed (see Fig. [5| . Similar ideas have been used recently in the design of sparse 
coding matrices for error-correcting codes [T5H2T]. 

Analysis of the performance of the seeded Belief Propagation procedure 

The s-BP procedure is based on the joint use of seeding measurement matrices and of the EM-BP message-passing 
reconstruction. We have studied it with two methods: direct numerical simulations and analysis of the performance 
in the large N limit. The analytical result was obtained by a combination of the replica method and of the cavity 
method (also known as 'density evolution' or 'state evolution'). The replica method is a standard method in statistical 
physics , which has been applied successfully to several problems of information theory [T71 EH EH] including 
compressed sensing [I2HI4] . It can be used to compute the free entropy function $ associated with the probability 
P(x) (see Appendix D), and the cavity method shows that the dynamics of the message-passing algorithm is a gradient 
dynamics leading to a maximum of this free-entropy. 

When applied to the usual case of the full F matrix with independent Gaussian-distributed elements (case L = 1), 
the replica computation shows that the free-entropy $(-D) for configurations constrained to be at a mean-squared 
distance D has a global maximum at D — when a > pq, which confirms that the Gauss-Bernoulli probabilistic 
reconstruction is in principle able to reach the optimal compression limit a = pQ. However, for aEM-BP > ck > Poj 
where aEM-BP is a threshold that depends on the signal and on the distribution P(x), a secondary local maximum of 
<I>(D) appears at D > (see Fig. [s]). In this case the EM-BP algorithm converges instead to this secondary maximum 
and does not reach exact reconstruction. The threshold q^em-bp is obtained analytically as the smallest value of a 
such that $(£)) is decreasing (Fig. [2|. This theoretical study has been confirmed by numerical measurements of the 
number of iterations needed for EM-BP to reach its fixed point (within a given accuracy). This convergence time of 
BP to the exact reconstruction of the signal diverges when a — > oem-bp (see Fig. |3|. For a < Oem-bp the EM-BP 
algorithm converges to a fixed point with strictly positive mean-squared error (MSE). This 'dynamical' transition is 
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FIG. 4: Construction of the measurement matrix F for seeded compressed sensing. The elements of the signal vector are split 
into L (here L = 8) equally-sized blocks, the number of measurements in each block is Mp = apN/L (here Qi = 1, Op = 0.5 
for p = 2, ... ,8). The matrix elements F^i are chosen as random Gaussian variables with variance Jq,p/N if variable i is in 
the block p and measurement fi in the block q. In the s-BP algorithm we use Jp,, = 0, except for Jp^p = 1, Jp,p-i ~ Ji, and 
Jp-i,p = J2. Good performance is typically obtained with relatively large Ji and small J2. 



similar to the one found in the cooling of Hquids which go into a super-cooled glassy state instead of crystallizing, and 
appears in the decoding of error correcting codes jl61 117) as well. 

We have applied the same technique to the case of seeding-measurement matrices {L > 1). The cavity method 
allows to analytically locate the dynamical phase transition of s-BP. In the limit of large N, the MSE Ep and the 
variance messages Vp in each block p ~ the density p, the mean 5;, and the variance of P(x) evolve 

according to a dynamical system which can be computed exactly (see Appendix E), and one can see numerically if 
this dynamical system converges to the fixed point corresponding to exact reconstruction (Ep = for all p). This 
study can be used to optimize the design of the seeding matrix F by choosing ai, L and Jp^q in such a way that the 
convergence to exact reconstruction is as fast as possible. In Fig. [3] we show the convergence time of s-BP predicted 
by the replica theory for different sets of parameters. For optimized values of the parameters, in the limit of a large 
number of blocks L, and large system sizes N/L, s-BP is capable of exact reconstruction close to the smallest possible 
number of measurements, a — po- In practice, finite size effects slightly degrade this asymptotic threshold saturation, 
but the s-BP algorithm nevertheless reconstructs signals at rates close to the optimal one regardless of the signal 
distribution, as illustrated in Fig. |2] 

We illustrate the evolution of the s-BP algorithm in Fig. [5] The non-zero signal elements are Gaussian with zero 
mean, unit variance, density po — 0-4 and measurement rate a = 0.5, which is deep in the glassy region where all 
other known algorithms fail. The s-BP algorithm first nucleates the native state in the first block and then propagates 
it through the system. We have also tested the s-BP algorithm on real images where the non-zero components of the 
signal are far from Gaussian, and the results are nevertheless very good, as shown in Fig. [T] This shows that the 
quality of the result is not due to a good guess of -P(x). It is also important to mention that the gain in performance 
in using seeding-measurement matrices is really specific to the probabilistic approach: we have computed the phase 
diagram of £1 minimization with these matrices and found that, in general, the performance is slightly degraded with 
respect to the one of the full measurement matrices, in the large N limit. This demonstrates that it is the combination 
of the probabilistic approach, the message-passing reconstruction with parameter learning, and the seeding design of 
the measurement matrix that is able to reach the best possible performance. 



Perspectives 

The seeded compressed sensing approach introduced here is versatile enough to allow for various extensions. One 
aspect worth mentioning is the possibility to write the EM-BP equations in terms of N messages instead of the M x N 
parameters as described in Appendix C. This is basically the step that goes from rBP 24J to AMP 11 algorithm. 
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FIG. 5: Evolution of the mean-squared error at different times as a function of the block index for the s-BP algorithm. Exact 
reconstruction first appears in the left block whose rate qi > oem-bp allows for seeded nucleation. It then propagates gradually 
block by block driven by the free entropy difference between the metastable and equilibrium state. After little more than 300 
iterations the whole signal of density po — 0.4 and size A'^ = 50000 is exactly reconstructed well inside the zone forbidden for 
BP (see Fig.[2|. Here we used L = 20, Ji = 20, J2 = 0.2, ai = 1.0, a = 0.5. 



It could be particularly useful when the measurement matrix has some special structure, so that the measurements 
y can be obtained in many fewer than M x N operations (typically in NlogN operations). We have also checked 
that the approach is robust to the introduction of a small amount of noise in the measurements (see Appendix H). 
Finally, let us mention that, in the case where a priori information on the signal is available, it can be incorporated 
in this approach through a better choice of (j), and considerably improve the performance of the algorithm. For signal 
with density pq, the worst case, that we addressed here, is when the non-zero components of the signal are drawn 
from a continuous distribution. Better performance can be obtained with our method if these non-zero components 
come from a discrete distribution and one uses this distribution in the choice of (j). Another interesting direction in 
which our formalism can be extended naturally is the use of non-linear measurements and different type of noises. 
Altogether, this approach turns out to be very efficient both for random and structured data, as illustrated in Fig. [T] 
and offers an interesting perspective for fpractical compressed sensing applications. Data and code are available online 

at [g. 
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Note: During the review process for our paper, we became aware of the work [39' in which the authors give a 
rigorous proof of our result, in the special case when pq — p and 4>q = 0, that the threshold a = po can be reached 
asymptotically by the s-BP procedure. 
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Appendix A: Proof of the optimality of the probabilistic approach 



Here we give the main lines of the proof that our probabihstic approach is asymptotically optimal. We consider the 
case where the signal s has iid components 



N 



Po{s) = - Po)S{si) + poMsi)] , (Al) 



with < po < 1- And we study the probability distribution 

N M / \ 

= ^ n ('^^^ - + '^'^(^')]) n E ^'''(^^ -'^n ' (^2) 

i=l ^=1 \ i / 

with a Gaussian 4>{x) of mean zero and unit variance. We stress here that we consider general 4>q, i.e. (f>Q is not 
necessarily equal to (^(x) (and po is not necessarily equal to p). The measurement matrix F is composed of iid 
elements F^j such that if /u belongs to block q and i belongs to block p then F^j is a random number generated from 
the Gaussian distribution with zero mean and variance Jg,p/N. The function Se{x) is a centered Gaussian distribution 
with variance e^. 

We show that, with probability going to one in the large N limit (at fixed a = M/N), the measure P (obtained 
with a generic seeding matrix F as described in the main text) is dominated by the signal ii a > po, a' > po (as long 
as (/'o(O) is finite). 

We introduce the constrained partition function: 

Y{D, e) = / n (dxi [(1 - p)S{xi) + p0(xi)]) [] ( E ^'^^(^^ " ) ^ ( E^^^ " *^)' > ' ^^^^ 

•' i=l M=l V i / \i=l / 

and the corresponding 'free entropy density': y{D,e) = limjv->.oo EpEg logy(D, e)/iV. The notations Ep and Eg 
denote respectively the expectation value with respect to F and to s. I denotes an indicator function, equal to one if 
its argument is true, and equal to zero otherwise. 

The proof of optimality is obtained by showing that, under the conditions above, lime_).o y{D, e)/[(a — po) log(l/e)] 
is finite if _D = (statement 1), and it vanishes if D > (statement 2). This proves that the measure P is dominated 
by D = 0, i.e. by the neighborhood of the signal Xi = si. The standard 'self-averageness' property, which states that 
the distribution (with respect to the choice of F and s) of logF(_D, e)/N concentrates around y{D, e) when N ^ oo, 
completes the proof. We give here the main lines of the first two steps of the proof. 

We first sketch the proof of statement 2. The fact that lim^^o y{D, e)/[{a — po) log(l/e)] = when D > can be 
derived by a first moment bound: 

y{D, e) < hm ^Eg log Y,^^{D, e) , (A4) 

iv— >cx> iV 

where Ya,nn{D,e) is the 'annealed partition function' defined as 

Y,UD,e)=E-pY{D,e). (A5) 

In order to evaluate Ya,nn{D,€) one can first compute the annealed partition function in which the distances between 
X and the signal are fixed in each block. More precisely, we define 

Z{n, • • • , ri, e) = Ef / n (dxi [(1 - p)5{xi) + pcl>{xi)]) JJ S, ( ^ F^i{xi - sM H M " ^ E " ^')' 

1=1 n=l \ i J p=l \ ieBp 

By noticing that the M random variables = J^i^^ii^i ~ ^i) are independent Gaussian random variables one 
obtains: 



L 



p=i \ 



9=1 



-Nap/2 



L 



JJ e^^'^)^^-"^ , (A6) 
p=i 



9 



where 



lim — log 

n— ^oo Ti 



/n / n 

n {dx^ [(1 - p)Six,) + d r - ^ ^(a 

4=1 \ 1 = 1 



(A7) 



The behaviour of ipir) is easily obtained by standard saddle point methods. In particular, when r — ?■ 0, one has 
tjj{r) ~ ipologr. 

Using (A6), we obtain, in the small e limit: 



lim -Eslogran„(^,e) 

A*— foo iV 



max 



1 1 

|^-^logr,--5: 



p=i 



L ^ 2 

p=i 



loe 



9=1 



(A8) 



where the maximum over ri, . . . , is to be taken under the constraint ri 



tl > LD. Taking the limit of e 



with a finite D, at least one of the distance must remain finite. It is then easy to show that 



lim -Eslogra„n(7^,e) 

N—>00 1\ 



log(l/e) 



a 



Po 



L 



{2a' - Po) 



(A9) 



where a' is the fraction of measurements in blocks 2 to L. As a' > poi this is less singular than log(l/e)(a — po); 
which proves statement 2. 

On the contrary, when _D = 0, we obtain from the same analysis 



lim 4Eslogra„n(i^,e) 



log(l/e)(a-po) 



(AlO) 



This annealed estimate actually gives the correct scaling at small e, as can be shown by the following lower bound. 
When 13 = 0, we define Vq as the subset of indices i where Si — 0, |Vo| = N{1 — pq), and Vi as the subset of indices i 
where Si ^ 0, |Vi| = Np^. We obtain a lower bound on Y{0, e) by substituting P(x) by the factors (1 — p)5{xi) when 
i (z Vq and p<j)[xi) when i G Vi. This gives: 



y (0, e) > exp{iV[(l - Po) log(l - p) + Po log(p) - {a/2) log(27r) + (po - a) log e)]} 



dui 



{si + tUi) exp 



ieVi 



(All) 



where Mij = X]J^=i FfiiFfij- The matrix Af , of size pqN x po-/V, is a Wishart-like random matrix. For a > po, generi- 
cally, its eigenvalues are strictly positive, as we show below. Using this property, one can show that, if Y[ 



) >0, 

the integral over the variables Ui in (All ) is strictly positive in the limit e 0. The divergence of 3^(0, e) in the limit 
e — >■ is due to the explicit term exp[Af(po — a)loge] in Eq. (All). 

The fact that all eigenvalues of M are strictly positive is well known in the case of L = 1 where the spectrum has 
been obtained by Marcenko and Pastur. In general, the fact that all the eigenvalues of M are strictly positive is 
equivalent to saying that all the lines of the aN x poA^ matrix F (which is the restriction of the measurement matrix 
to columns with non-zero signal components) are linearly independent. In the case of seeding matrices with general 
L, this statement is basically obvious by construction of the matrices, in the regime where in each block q, Uq > po 
and Jqq > 0. A more formal proof can be obtained as follows. We consider the Gaussian integral 



/n 
Y\ d(f>i exp 



2 ^ y ^ 



(A12) 



This quantity is finite if and only if v is smaller than the smallest eigenvalue, Amin, of M. We now compute the 
annealed average Zs^nn{v) = EfZ{v). If Zs^nn{v) is finite, then the probability that Amin < v goes to zero in the large 
N limit. Using methods similar to the one above, one can show that 



2L 



logEF^(w) = max 

n ri,...,rp 



.p=i 



vrp) 



\ q 



(A13) 
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The saddle point equations 



1 ^1 \ ^ Uq Jqp 



(AM) 



qs' s 



have a solution Sit v — (and by continuity also at t' > small enoug h), when i J2q=i ^ = > Mit 

can be found 



Po 



for instance by iteration). Therefore, ^ logEF^{w) is finite for some v > small enough, and therefore Amin > 0. 



Appendix B: Derivation of Expectation maximization Belief Propagation 



In this and the next sections we present the message-passing algorithm that we used for reconstruction in compressed 
sensing. In this section we derive its message-passing form, where 0{NM) messages are being sent between each 
signal component i and each measurement /i. This algorithm was used in |24j . where it was called the relaxed belief 
propagation, as an approximate algorithm for the case of a sparse measurement matrix F. In the case that we use 
here of a measurement matrix which is not sparse (a finite fraction of the elements of F is non-zero, and all the 
non-zero elements scale as 1/y/N), the algorithm is asymptotically exact. We show here for completeness how to 
derive it. In the next section we then derive asymptotically equivalent equations that depend only on 0{N) messages. 
In statistical physics terms, this corresponds to the TAP equations [H] with the Onsager reaction term, that are 
asymptotically equivalent to the BP on fully connected models. In the context of compressed sensing this form of 
equations has been used previously [11^ and it is called approximate message passing (AMP). In cases when the matrix 
F can be computed recursively (e.g. via fast Fourier transform), the running time of the AMP-type message passing 
is 0{N log N) (compared to the 0{NM) for the non-AMP form). Apart for this speed-up, both classes of message 
passing give the same performance. 

We derive here the message-passing algorithm in the case where measurements have additive Gaussian noise, the 
noiseless case limit is easily obtained in the end. The posterior probability of x after the measurement of y is given 
by 



P(x) 



1 ^ 



M 



1 



(Bl) 



where Z is a normalization constant (the partition function) and is the variance of the noise in measurement /z. 
The noiseless case is recovered in the limit — ^ 0. The optimal estimate, that minimizes the MSB with respect to the 
original signal s, is obtained from averages of Xi with respect to the probability measure ^'(x). Exact computation of 
these averages would require exponential time, belief propagation provides a standard approximation. The canonical 
BP equations for probability measure P(x) read 



ix.j)dxj 



^ [(1 - p)S{x,) + p(l){xi)] Y[ m^^i{x,) . 



(B2) 
(B3) 



where Z^^^ and Z^^^ are normalization factors ensuring that J da;iTO^_>i(xj) = J dximi^^{xi) = 1. These are 
integral equations for probability distributions that are still practically intractable in this form. We can, however, 
take advantage of the fact that after proper rescaling the linear system y = Fx is such a way that elements of y 
and x are of 0(1), the matrix has random elements with variance of 0(1/N). Using the Hubbard-Stratonovich 
transformation 



1 



dAe 2A+ A 



(B4) 



for Lo = i^j^iF^jXj) we can simplify eq. (B2) as 

1 



dAe 



n 



axjmj^fj_[Xj)e ' ^' * 



(B5) 
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Now we expand the last exponential around zero because the term _F)jj is small in N, we keep all terms that are of 
0{1/N). Introducing means and variances as new messages 



dx^ X.^ TTii^^i^Xi) 



we obtain 

1Tlfj,—^i{Xi) 



1 



Performing the Gaussian integral over A we obtain 

1 



e Y[ 



^fi^i {Xi ) 



2tt 



A 



where we introduced 



B 



(B6) 
(B7) 

(B8) 

(B9) 

(BIO) 
(Bll) 



and the normalization Z^^^ contains all the Xi-independent factors. The noiseless case corresponds to = 0. 
To close the equations on messages a^-j.^ and Vi^^ we notice that 

m,^^{x,) = -J— [(1 - p)S{x,) + pcj,(x,)] e-'^^-<*^ A-,^,+x, E,^, b,^, ^3^2) 

Messages a.i^^ and w^-i.^ are respectively the mean and variance of the probability distribution mi^fj_(xi). For general 
(j){xi) the mean and variance ( B6||B7 ) will be computed using numerical integration over Xi. Eqs. ( B6|B7 ) together 
with ( B10|B11 1 and (B12) then lead to closed iterative message-passing equations. 

In all the specific examples shown here and in the main part of the paper we used a Gaussian 4>{xi) with mean x 
and variance cr^. We define two functions 



fa{X,Y) = 

fc{X,Y) = 



p{Y + x/a^) 
a(l/a2 +X)3/2 

P 



(1 - p)e 2(l/<,2+x)^2<,2 ^ 



(B13) 



1 



_a(l/(72 +X)-V2 V 1/(72 + X 
Then the closed form of the BP update is 



[1 — p)e 2(1/,t2+x) "^2^2 



(7(1/(72+^)1/2 



fa{X,Y). 



77^ A" 



^ J 0,1 — fa ^ ^7— »ii ^ ^ ^7— >-t^ ) 

= /c ( ^ ^ ^7-i.ij ^ ^ B~f^i j , 
\ 7 7 / 



(B14) 

(B15) 
(B16) 



where the and are the mean and variance of the marginal probabilities of variable Xi . 

As we discussed in the main text the parameters p, x and a are usually not known in advance. However, their values 
can be learned within the probabilistic approach. A standard way to do so is called expectation maximization [34) . 
One realizes that the partition function 



z{p,x,g) = / n^^«n 



AT 



N r 



i=l 1=1 



{l-p)5{x,) + -^e 2.2 



M ^ 



■3 2Au Si=l ^[ii'^i 



(B17) 
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is proportional to the probability of the true parameters po,s,cro, given the measurement y. Hence to compute the 
most probable values of parameters one searches for the maximum of this partition function. Within the BP approach 
the logarithm of the partition function is the Bethe free entropy expressed as[T7] 



F{p, x,a)^Y. + E - E log 



where 



dj:j]^mp^,j(.Xj) 

rid- n 



{1 - p)6{x^) + - 



[xi)- 



/27rCT 



<:\ximf^^i{xi)'mi^fj_{xi) . 



The stationarity conditions of Bethe free entropy (B18) with respect to p leads to 



P ■ 



(Vj+x/o-01_^ 
) 2(l/cr^ + (7j) 2tT 



cr(l/o-2+C/,)2 

where C/^ — ^7->.i, and = B^^i. Stationarity with respect to x and cr gives 

^ _ E»a* 



(Vi+x/,T2)2 _2 



p+ (1 - p)cr(l/cr2 + (7^)2 6 2<i/-2+^.) 



/9+(l-p)(T(l/cr2 + J7,)5e ^d/" 



(Vi+!ir/<T2)2 _2 



(B18) 

(B19) 

(B20) 
(B21) 

(B22) 



(B23) 



(B24) 



In statistical physics conditions ( B22 1 are known under the name Nishimori conditions [351 ISZ] • In the expectation 
maximization eqs. ( B22|B24 | they are used iteratively for the update of the current guess of parameters. A reasonable 
initial guess is pinit. — ol. The value of p^s can also be obtained with a special line of measurement consisting of a 
unit vector, hence we assume that given estimate of p the x = pa's/ p. In the case where the matrix F is random 
with Gaussian elements of zero mean and variance we can also use for learning the variance: '^^^iVfi/N — 



Appendix C: AMP-form of the message passing 

In the large N limit, the messages ai_>^ and Vi^^ are nearly independent of /i, but one must be careful to keep the 
correcting Onsager reaction terms. Let us define 

i i 



Then we have 



V. = E + ^ E + (u., E p% -A- (C4) 

y + 7^ - F^^^v,^^ ^ A^ + ^ A^ + 7,, 
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We now compute uj^ 



fa {Ui - A^^i,Vi - Bf_,^i) ~ Ui 



(C5) 



To express oj^ = J2i ^t^i'^i-^fj.' '^^ that the first correction term has a contribution in F^^, and can be safely 
neglected. On the contrary, the second term has a contribution in F^^^ which one should keep. Therefore 



-7m 



The computation of 7^ is similar, it gives: 



7m = V Pu^V^ 



> F-" 



dfc 



A, 



-7m 



dY 



\F%UU,,Vi). 



(C6) 



(C7) 



For a known form of matrix F these equations can be slightly simplified further by using the assumptions of the 
BP approach about independence of F^i and BP messages. This plus a law of large number implies that for matrix F 



with Gaussian entries of zero mean and unit variance one can effectively 'replace' every F'^^ by 1/A^ in eqs. (C3 C4) 



and ( C6|C7|). T his leads, for homogeneous or bloc matrices, to even simpler equations and a slightly faster algorithm 



1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 



^ 7t 

Eqs. ( C3|C4 ) and ( C6|C7 ), with or without the later simplification, give a system of closed equations. They are a 
special form (P(x) and hence functions /a, /c are different in our case) of the approximate message passing of |11) . 

The final reconstruction algorithm for general measurement matrix and with learning of the P(x) parameters can 
hence be summarized in a schematic way: 

EM-BP(y^, F^i, criterium, ^max) 

Initialize randomly messages Ui from interval [0, 1] for every component; 
Initialize randomly messages Vi from interval [—1, 1] for every component; 
Initialize messages uip_ y^; 

Initialize randomly messages 7^ from interval ]0, 1] for every measurement; 
Initialize the parameters p a, x 0, cr^ -s— 1. 
conv criterium +1; i 0; 
while conv > criterium and t < imax^ 
dot^t+1; 

for each component i: 



do X°'d^/a(i7„l/,); 

Update Ui according to eq. 

Update Vi according to eq. 
for each measurement fi: 

do Update according to eq. 

Update 7^ accordi ng to eq. 
Update p according to eq. (B22). 



Update X and according to eq. ( |B24[ ). 
for each component i: 

do ^ fa{Ui, Vi); 
conv ^ mean(|a;i — x"'*^!); 
return signal components x 

Note that in practice we use 'damping' (at each update, the new message is obtained as u times the old value plus 
1 — u times the newly computed value, with a damping < u < 1, typically u = 0.5) for both the update of messages 
and learning of parameters, empirically this speeds up the convergence. Note also that the algorithm is relatively 
robust with respect to the initialization of messages. The reported initialization was used to obtain the results in 
Fig. 1 of the main text. However, other initializations are possible. Note also that for specific classes of signals or 
measurement matrices the initial conditions may be adjusted to take into account the magnitude of the values F^i 
and y^. 

For a general matrix F one iteration takes 0{NM) steps, we observed the number of iterations needed for conver- 
gence to be basically independent of N , however, the constant depends on the parameters and the signal, see Fig. 3 
in the main paper. For matrices that can be computed recursively (i.e. without storing all their NM elements) a 
speed-up is possible, as the message-passing loop takes only 0{M -t- N) steps. 
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Appendix D: Replica analysis and density evolution: full measurement matrix 

Averaging over disorder leads to replica equations that are describing the — > oo behavior of the partition function 
as well as the density evolution of the belief propagation algorithm. The replica trick evaluates EF,s(log^) via 

$ = -E(log Z)^- hm . (Dl) 

jy jy n~>Q n 

In the case where the matrix F is the full measurement with all elements independent identically distributed from 
a normal distribution with zero mean and variance unity, one finds that $ is obtained as the saddle point value of the 
function: 

^Q,q,m,Q,q,m) = --^ A + Q-q -log (A + Q - q) + ^ - mm + ^ 

+ J'DzJ ds[(l-po)'^(s) + Po0o(s)]log|y d:Ee-^-'+™-^+^^-[(l-p),5(x)+p0(x)]|. (D2) 

Here Vz is a Gaussian integration measure with zero mean and variance equal to one, po is the density of the signal, 
and 0o(s) is the distribution of the signal components and (s^) — J dss'^(f>Q{s) is its second moment. A is the variance 
of the measurement noise, the noiseless case is recovered by using A = 0. 
The physical meaning of the order parameters is 

Q'^j^^i^t)^ ^""^X!^^*^^' ^=j^^s^i^i)- (D3) 

i i i 

Whereas the other three m, q, Q are auxiliary parameters. Performing saddle point derivative with respect to 
m,q,Q — q,m,q,Q + q we obtain the following six self-consistent equations (using the Gaussian form of 4>{x), with 
mean x and variance a^): 

- " -^'-Z- ^ g-2m + p,(,s^)+A ^ ^^^^ 



A + Q-q ^ ^' ^ A + Q-q (A + Q-q) 



P°P f T, f A A t \ TTis + zy/q + x/a^ 

2(Q + 9 + 1/<t2) I „ 



Q + q + l/a^J J {l-p)aJQ + q+l/a^e 



(l-Po)p z^/q + x/a'^ 

' Vz z 



{Q + q + l/(7^)VqJ (1- p)a^Q + q + l/cr^e^^'' 2lQlT+i%'^) + p 



Pop f A A f \ ms + z^ + x/a^ 



{Q + q + l/(7^)VqJ J {1- p)a^Q + q + l/cr^eJ''' 'TtO+Z+v^'^)' +p 



i'^-Po)P f^^_ [zy/q + x/a^f + Q + q + l/a^ 



[LJ + q + i/a) J {I- p)cr^Q + q + l/a^e^ 2(q+«+i/.2, _^ ^ 

Pop [ T, [ A A ( \ {rhs + z^ + x/a'^f + Q + q+\/a'^ 

/ Vz j ds0o(s) — 5 — , - ^ ,,2 . (D7) 

{Q + q + l/a-^yJ J ■ ^ r. ^ — ^_ (^.+.v^+x/^ ->^ 



[l- p)a\J Q + q + l/a'^e'^ 2(<5+4+i/<.^) + p 

We now show the connection between this replica computation and the evolution of belief propagation messages, 
studying first the case where one does not change the parameters p, x and a. Let us introduce parameters mep, gep, 
Qbp defined via the belief propagation messages as: 



N N N 

i—l i—1 i—1 



The density (state) evolution equations for these parameters can be derived in the same way as in |11|, I33j , and this 
leads to the result that ttt-bp , (/bp j Qbp evolve under the update of BP in exactly the same way as according to iterations 
of eqs. ( D5|D7 ). Hence the analytical eqs. ( D5|D7 ) allow to study the performance of the BP algorithm. Note also 
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that the density evolution equations are the same for the message-passing and for the AMP equations. It turns out 
that the above equations close in terms of t wo para meters, the mean-squared error E'gp = g^p — 2r7igp -I- po{s'^) and 



the variance V, 



(*) 



BP 



?BP - 9bp- From eqs. (|D4 



D7) easily gets a closed mapping (^E^p^\v^p^^^ = / (^^bP' ^b*p) ■ 
In the main text we defined the function $(-0) which is the free entropy restricted to configurations x for which 
D = J2i=ii^i ~ Si)'^/N is fixed. This is evaluated as the saddle point over Q, q, Q, q, rh of the function ^{Q, q, {Q ~ 
D + /Oo(s^))/2, Q, q, rh). This function is plotted in Fig. 3(a) of the main text. 

In presence of Expectation Maximization learning of the parameters, the density evolution for the conditions (B22 ) 
and (|B24|) are 



p{t+i) ^ pit) 



Vz / da;o [(1 - Po)5{xq) + po0o(a;o)] 



2N(t+l) ^ 



g{Q + q,mxo + z^/q) \ 
I- p + pg{Q + q, mxo + zy^) J 

^ y 

1 - P + PaiQ + q, nixa + z^/q) ) 
j Vz J dxo [(1 - po)5{xo) + pocjJoixo)] fa{Q + q,mxo + z^q) 

T, f A ^^^ \xt \ , M giQ + q, 171X0 + z^/q) \ 

Vz / dxo[{l - po)o(xo) + poM^o)] 77; — 7=7 

J 1 - p + pg[Q + q, niXQ + z^q) j 

I 'Dz J dxo [(1 - Pq)S{xo) + PoMxo)] [fa{Q + q,mxo + z^qf + /^(Q + q^rnx^ + z^q)\ 

g{Q + q.mxa + z^fq) 



(D9) 



(DIO) 



Vz \ Axo{(\ - pa)^(xQ) + pq(1)q(xo)\ 



i;(*+l)l2 



1 - p + pg{Q + q, mxo + 



The density evolution equations now provide a mapping 



J^{t+l) y{t + l} (t + 1) (f + 1) (t + l) 

, EM-BP' EM-BP ' '-^ 



At+i) 



J \ -^EM-BP' '^EM-BP' '-^ 
it) 



obtained by complementing the previous equations on 



it) 



(Dll) 



(D12) 



EM-BP' EM-BP 



with the update equations (D9 



DIO 



Dill. 



The next section gives explicitly the full set of equations in the case of seeding matrices, the ones for the tull matrices 
are obtained by taking L = 1. These are the equations that we study to describe analytically the evolution of EM-BP 
algorithm and obtain the phase diagram for the reconstruction (see Fig. 2 in the main text). 



Appendix E: Replica analysis and density evolution for seeding-measurement matrices 

Many choices of Ji and J2 actually work very well, and good performance for seeding-measurement matrices can be 
easily obtained. In fact, the form of the matrix that we have used is by no means the only one that can produce the 
seeding mechanism, and we expect that better choices, in terms of convergence time, finite-size effects and sensibility 
to noise, could be unveiled in the near future. 

With the matrix presented in this work, and in order to obtain the best performance (in terms of phase transition 
limit and of speed of convergence) one needs to optimize the value of Ji and J2 depending on the type of signal. 
Fortunately, this can be analysed with the replica method. The analytic study in the case of seeding measurement 
matrices is in fact done using the same techniques as for the full matrix. The order parameters are now the MSE 
Ep = qp — 2mp + po{s'^) and variance = Qp — qp in each block p G {1, ... ,1/}. Consequently, we obtain the final 
dynamical system of 2L -I- 3 order parameters describing the density evolution of the s-BP algorithm. The order 
parameters at iteration t + 1 of the message-passing algorithm are given by: 
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E;(t+i) 
q 



y(t+i) 



da;0 [(1 - po)'5(xO) + po0o(a^")] j Vz (^fa (q, + Qg, m,xO + 
dx° [(1 - po)'5(a;") + paMx°)] / 2?^/c (q? + Qq, rhqX° + z^?^) 



(El) 
(E2) 



P + P.9(Qp + ^p, rhpXa + ^:v^) 



/Pz I dxQ [(1 - po)'5(a;o) + /3o'/>o(a;o)] - 
V p=i 1 

1 1 dxo [(1 - po)^(a;o) +Po0o(a;o)] :^ 77^ — : 7=r 

^ p=iJ J 1- p + pg(Qp + qp,mpXo + Zy^qp) 



(E3) 



p u 



p=l 



^ jDz dxoiil - Po)Sixo)+ po<j)oixo)]fa{Qp + qp,rhpXo + zy%) 



[ T^z j dxo [{1 - po)S{xo) + po(l>o{xo)] — 



g{Qp + qp, rhpXQ + z^p) 



p \L 



E 

p=l 



- p + pgiQp + qp,rhpXo + Zy/qp) J 

Vz j dxo [(1 - Po)(5(xo) + Po0o(a;o)] [/a(Qp + qp^rhpXo + zy%)'^ + fdQp + qp,rhpXQ + z^/Yp) 

, -1 

g{Qp + gp, mpXo + z^p) \ 



(E4) 



jS^ j T^z j dxo [[I - po)5[xf)) + pf)(l)o{xo)] 
p=i •' •' 



F;(*+1)i2 



(E5) 



where: 



99 



^E 



JpqOlp 



(t) 



- _ 1 \ - a„J, 

- tZ^ 



J^p'Jpq 



p " 



99 = 



99- 



The functions /a(^, l'), fciX,Y) were defined in ( B13||B14 ), and the function g is defined as 



=exp 



(Y + x/f72)2 



(E6) 
(E7) 
(E8) 

(E9) 



This is the dynamical system that we use in the paper in the noiseless case (A = 0) in order to optimize the values 
of ai, Ji and J2 . We can estimate the convergence time of the algorithm as the number of iterations needed in 
order to reach the successful fixed point (where all Ep and Vp vanish within some given accuracy). Figure [6] shows 
the convergence time of the algorithm as a function of Ji and J2 for Gauss-Bernoulli signals. 

The numerical iteration of this dynamical system is fast. It allows to obtain the theoretical performance that can be 
achieved in an infinite- A'^ system. We have used it in particular to estimate the values of L, ai, Ji, J2 that have good 
performance. For Gauss-Bernoulli signals, using optimal choices of Ji, J27 we have found that perfect reconstruction 
can be obtained down to the theoretical limit a = p^hy taking L — > 00 (with correction that scale as 1/i). Recent 
rigorous work by Donoho, Javanmard and Montanari j?^ extends our work and proves our claim that s-BP can reach 
the optimal threshold asymptotically. 

Practical numerical implementation of s-BP matches this theoretical performance only when the size of every block 
is large enough (few hundreds of variables). In practice, for finite size of the signal, if we want to keep the block-size 
reasonable we are hence limited to values of L of several dozens. Hence in practice we do not quite saturate the 
threshold a = po, but exact reconstruction is possible very close to it, as illustrated in Fig. 2 in the main text, where 
the values that we used for the coupling parameters are listed in Table |lj 

In Fig. 3, we also presented the result of the s-BP reconstruction of the Gaussian signal of density po = 0.4 for 
different values of L. We have observed empirically that the result is rather robust to choices of Ji, J2 and ai. 
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Ji 

FIG. 6: (color online): Convergence time of the s-BP algorithm with L = 2 as a function of Ji and J2 (in log scale) for 
Gauss-Bernoulli signal with po — 0.1. The color represents the number of iterations such that the MSB is smaller than 10~*. 
The white color region gives the fastest convergence. The measurement density is fixed to a = 0.25. The parameters in s-BP 
are chosen as L = 2 and ai — 0.3. The axes are in logj^g scale. 
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0.16 


20 




0.4 


0.441 


0.8 


0.422 


16 


0.16 


20 


0.6 


0.624 


0.9 


0.609 


4 


0.04 


20 




0.6 


0.630 


0.95 


0.613 


16 


0.16 


20 


0.8 


0.816 


0.95 


0.809 


4 


0.04 


20 




0.8 


0.820 


1 


0.811 


4 


0.04 


20 



TABLE I: Parameters used for the s-BP reconstruction of the Gaussian signal (left) and of the binary signal (right). 

In this case, in order to demonstrate that very different choices give seeding matrices which are efficient, we used 
ai = 0.7, and then Ji = 1043 and J2 = lO"'' for L = 2,Ji = 631 and J2 = 0.1 for L = 5, Ji = 158 and J2 = 4 for 
L = 10 and Ji — 1000 and J2 = 1 for L = 20. One sees that a common aspect to all these choices is a large ratio 
J1/J2. Empirically, this seems to be important in order to ensure a short convergence time. A more detailed study of 
convergence time of the dynamical system will be necessary in order to give some more systematic rules for choosing 
the couplings. This is left for future work. 

Even though our theoretical study of the seeded BP was performed on an example of a specific signal distribution, 
the examples presented in Figs. 1 and 2 show that the performance of the algorithm is robust and also applies to 
images which are not drawn from that signal distribution. 

Appendix F: Phase diagram in the variables used by fE] 

We show in Fig. [t] the phase diagram in the convention used by [5] , which might be more convenient for some 
readers. 
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FIG. 7: (color online): Same data as in Fig. 2 in the main paper, but using the convention of [5]. The phase diagrams is now 
plot as a function of the under-sampling ratio pnr = K/M — p/a and of the over-sampling ratio Sdt = M/N = a. 



Appendix G: Details on the phantom and Lena examples 

In this section, we give a detailed description of the way we have produced the two examples of reconstruction in 
Fig. 1. It is important to stress that this figure is intended to be an illustration of the s-BP reconstruction algorithm. 
As such, we have used elementary protocols to produce true iiT-sparse signals and have not tried to optimize the 
sparsity nor to use the best possible compression algorithm; instead, we have limited ourselves to the simplest Haar 
wavelet transform, to make the exact reconstruction and the comparison between the different approaches more 
transparent. 

The Shepp-Logan example is a 128^ picture that has been generated using the Matlab implementation. The Lena 
picture is a 128^ crop of the 512^ gray version of the standard test image. In the first case, we have worked with 
the sparse one-step Haar transform of the picture, while in the second one, we have worked with a modified picture 
where we have kept the 24 percent of largest (in absolute value) coefficients of the two-step Haar transform, while 
putting all others to zero. The datasets of the two images are available online [7 . Compressed sensing here is done 
as follows: The original image is a vector o of = pixels. The unknown vector x — Wo are the projections of 
the original image on a basis of one- or two-steps Haar wavelets. It is sparse by construction. We generate a matrix 
F as described above, and construct G = FW. The measurements are obtained by y = Go, and the linear system 
for which one does reconstruction is y = Fx. Once x has been found, the original image is obtained from o ~ W^^x. 
We used EM-BP and s-BP with a Gauss-Bernoulh P(x). 

On the algorithmic side, the EM-BP and £i experiments were run with the same full gaussian random measurement 
matrices. The minimization of the £i norm was done using the £i-magic tool for Matlab, which is a free implementation 
that can be download at http:/ /users. ece.gatech.edu/~ just in/ 11 magic/, using the lowest possible tolerance such 
that the algorithm outputs a solution. The coupling parameters of s-BP are given in Table |IT] A small damping was 
used in each iterations, we mixed the old and new messages, keeping 20% of the messages from the previous iteration. 
Moreover, since for both Lena and the phantom, the components of the signal are correlated, we have permuted 
randomly the columns of the sensing matrix F. This allows to avoid the (dangerous) situation where a given block 
contains only zero signal components. 

Note that the number of iterations needed by the s-BP procedure to find the solution is moderate. The s-BP 
algorithm coded in Matlab finds the image in few seconds. For instance, the Lena picture at a = 0.5 requires about 
500 iterations and about 30 seconds on a standard laptop. Even for the most difficult case of Lena at a = 0.3, we need 
around 2000 iterations, which took only about 2 minutes. In fact, s-BP (coded in Matlab) is much faster than the 
Matlab implementation of £i-magic on the same machine. We report the MSE for all these reconstruction protocols 
in Table Uni 
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a 




a 


Ji 


J2 


L 




a 


ai 


a' 


Ji 


J2 


L 


0.5 


0.6 


0.495 


20 


0.2 


45 




0.6 


0.8 


0.585 


20 


0.1 


30 


0.4 


0.6 


0.395 


20 


0.2 


45 




0.5 


0.8 


0.485 


20 


0.1 


30 


0.3 


0.6 


0.295 


20 


0.2 


45 




0.4 


0.8 


0.385 


20 


0.1 


30 


0.2 


0.6 


0.195 


20 


0.2 


45 




0.3 


0.8 


0.285 


20 


0.1 


30 


0.1 


0.3 


0.095 


1 


1 


30 




0.2 


0.5 


0.195 


1 


1 


30 



TABLE II: Parameters used for the s-BP reconstruction of the Shepp-Logan phantom image (left) and of the Lena image 
(right). 





a = 0.5 


a = 0.4 


a = 0.3 


Q = 0.2 


a = 0.1 







0.0055 


0.0189 


0.0315 


0.0537 


BP 








0.0213 


0.025 


0.0489 


s-BP 














0.0412 





a = 0.6 


a = 0.5 


a = 0.4 


a = 0.3 


a = 0.2 







4.48.10"'' 


0.0015 


0.0059 


0.0928 


BP 





4.94.10"'' 


0.0014 


0.0024 


0.0038 


s-BP 














0.0038 



TABLE III: Mean-squared error obtained after reconstruction for the Shepp-Logan phantom (left) and for Lena (right) with 
£i minimization, BP and s-BP. 
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Noisy CS with Gauss-Bernoulli (Pq=0.2) signals 
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FIG. 8: (color online): Mean-squared error as a function of a for different values of the measurement noise. Left: Gauss- 
Bernoulli signal with po = 0.2, and measurement noise with standard deviation y/~A = 10"'*. The EM-BP (L — 1) and s-BP 
strategy {L = 4, Ji = 20, J = 0.1, ai — 0.4) are able to perform a very good reconstruction up to much lower value of a than 
the £i procedure. Below a critical value of a, these algorithms show a first order phase transition to a regime with much larger 
MSE. Right: Gauss-Bernoulli signal with po ~ 0.4, with a noise with standard deviation \/A = 10""*, 10"*, 10"^. s-BP, with 
L = 9, Ji — 30, J2 = 8, ai = 0.8 decodes very well for all these values of noises. In this case, £1 is unable to reconstruct for all 
a < 0.75, well outside the range of this plot. 



Appendix H: Performance of the algorithm in the presence of measurement noise 

A systematic study of our algorithm for the case of noisy measurements can be performed using the replica analysis, 
but goes beyond the scope of the present work. In this section we however want to point out two important facts: 
1) the modification of our algorithm to take into account the noise is straightforward, 2) the results that we have 
obtained are robust to the presence of a small amount of noise. As shown in section |B] the probabiUty of x after the 
measurement of y is given by 

i=l fi=l V M 

is the variance of the Gaussian noise in measurement /i. For simplicity, we consider that the noise is homogeneous, 
i.e., = A, for all /x, and discuss the result in unit of standard deviation \/A. The AMP-form of the message 
passing including noise have already been given in section [C] The variance of the noise, A, can also be learned via 
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expectation maximization approach, which reads: 

A ^ ^ (H2) 

The dynamical system describing the density evolution has been written in Sect.|Ej It just needs to be complemented 
by the iteration on A according to Eq. ( |H2[ ) . We have used it to study the evolution of MSE in the case where both 
P{x) and signal distribution are Gauss-Bernoulli with density p, zero mean and unit variance for the full measurement 
matrix and for the seeding one. Figure [8] shows the MSE as a function of a for a given A, using our algorithm, and the 
£i minimization for comparison. The analytical results obtained by the study of the dynamical system are compared 
to numerical simulations, and agree very well. 
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